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1.  INTRODUCTION 

Much  time  and  effort  have  been  expended  in  attempting  to  approximate  the  Fourier  transform  of 
functions.  The  Fourier  transform,  h,  of  a function  g is  given  by 

h(f)=[  g(t)ep27riftdt,  p = ± i.  (1) 

In  addition,  for  reasonable  functions,  the  inverse  to  this  transform  is  given  by 

g(t)=jf  h(0e“p27riftdf.  (2) 

For  this  report,  let  g be  real  and  continuous  and  vanish  outside  the  interval  [0,T] . In  this  case,  h(f)  is 
continuous,  and  g is  the  inverse  transform  of  h.  The  definition  of  the  Fourier  transform  may  easily  be  ex- 
tended to  distributions  of  which  the  Dirac  6 function  is  an  example. 

This  report  discusses  four  different  transforms  that  under  certain  conditions  may  be  considered  ap- 
proximations to  the  Fourier  transform.  The  Discrete  Fourier  Transform  (DFT)  and  the  Fast  Fourier 
Transform  (FFT)  are  well  known  and  represent  the  transform  of  impulses;  the  other  two  are  the  transforms 
of  piece-wise  linear  functions  and  are  implemented  by  the  subprograms  FLAT  and  NUFT. 

The  DFT  replaces  the  original  function,  g,  with  N impulses  (point  masses  at  N equispaced  points  in 
[0,  T]  weighted  by  the  value  of  g at  these  points).  The  Fourier  transform  of  the  sum  of  these  impulses  is 
then  computed.  Thus,  the  DFT  has  a value  at  each  frequency  point.  The  FFT  produces  a sampling  of  the 
DFT  at  N equispaced  frequency  points.  The  FFT  makes  use  of  the  Cooley-Tukey  algorithm  to  calculate 
these  values  rapidly  and  efficiently. 

Subroutine  NUFT  is  an  algorithm  that  calculates  the  Fourier  transform  of  virtually  any  piece-wise 
linear  function.  If  the  endpoints  of  these  linear  segments  are  equispaced,  N (tire  number  of  points)  is  a 
power  of  2,  and  tire  function  vanishes  outside  the  interval,  subroutine  FLAT  may  be  used  to  sample  NUFT 
at  N equispaced  frequency  points.  Subroutine  FLAT  makes  use  of  the  Cooley-Tukey  algorithm  to  produce 
these  results  in  about  the  same  computer  time  as  an  FFT. 

Since  each  of  these  transforms  is  the  Fourier  transform  of  an  approximation  to  the  original  function, 
the  degree  to  which  each  agrees  witli  the  Fourier  transform  is  largely  determined  by  the  goodness  of  fit  to 
tiic  original  functions. 

Since  it  is  of  great  interest  to  apply  these  transforms  to  the  analysis  of  rapidly  varying  transient 
signals,  this  report  discusses  also  the  routines  used  to  process  these  signals  at  the  Harry  Diamond  Labora- 
tories (HDL).  Also  discussed  is  Mimipulse  a function  used  to  simulate  the  response  of  systems  to  electro- 
magnetic pulse  (liMP) -as  well  as  the  calculation  of  its  transform.  Mimipulse  was  often  employed  as  a means 
of  testing  the  validity  of  the  various  processes. 

The  programs  are  listed  in  appendix  A. 


Preceding  page  blank 
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2.  THE  DISCRETE  FOURIER  TRANSFORM 


Let  At  > 0 and  N be  a positive  integer;  then  the  DFT  of  g of  order  N and  increment  At  is  given  by 

h(0=NI  g(jAt)ep2*ifAt,p  = ±l  . (3) 

j=0 

This  transform  is  of  interest  because  it  and  its  inverse  are  relatively  easy  to  compute. 

Insight  into  the  DFT  may  be  gained  by  observing  that  if  we  define 


N-l 


gj(tJ=  1 g(jAt)5(t  - jAt) , 
j=0  . 


(4) 


I N-l 

where  6 = Dirac  delta  function.  That  is,  gj  is  a sum  of  impulses  at  the  points  jjAt  , and  the  DFT  of  g 
is  identical  to  the  Fourier  transform  of  gj . Also,  the  DFT  defines  a mapping  from  functions  defined  on  a 
discrete  set, 


N-l 


I j=0 


to  functions  defined  for  all  f,  -«  < f < «,  for  if  two  functions  agree  on  this  discrete  set,  their  DFT’s  agree 
everywhere. 

3.  THE  FAST  FOURIER  TRANSFORM 


The  FFT  for  certain  N is  an  extremely  fast  method  of  sampling  the  DFT.  The  heart  of  the  F'F'T  is  the 
Cooley-Tukey  algorithm,  which  is  a very  efficient  means  of  summing  the  scries 


N-l 

>. 

j«0 


p2»«k/N 

Jj° 


,P°  *1 


(5) 


for  integer  k,  when  N is  a composite  number.  The  (,'ooloyTukcy  algorithm  lias  highest  efficiency  when 
N = 2m.  In  that  case,  the  ratio  ol  computer  time  spent  summing  by  use  of  this  algorithm,  as  opposed  to 
more  conventional  means,  is  approximately  m/N.  The  Cooley-Tukey  algorithm  is  applied  to  the  DFT  by  the 
observation  that  equation  (3)  reduces  to 


h(kAI)  ° 


N_| 

/. 


gOAth'1’ 


2 rrijk/ N 


(<>) 


if  f is  replaced  by  kAf,  where  Al  “ l/(NAt)  . 
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The  output  from  the  FFT  is  an  ordered  array  of  N complex  numbers.  The  first  (N/2  + 1)  numbers  of 
the  array  represent  the  values  at  frequencies  in  [0,  (N/2)Af)  inclusive  while  the  rest  of  the  array  represents 
the  values  at  frequencies  in  [(-N/2  + l)Af,  -Af]  inclusive.  Thus,  in  this  array,  the  negative  frequency  values 
follow  the  positive  ones,  and  the  highest  frequency  represented  is  (N/2)Af.  For  a real  function  g,  the  value 
at  a negative  frequency  is  the  complex  conjugate  of  the  value  at  the  corresponding  positive  frequency. 

3.1  Aliasing-Theorem  1 f Cooley) 

The  key  to  understanding  aliasing  is  a remarkable  theorem  of  Cooley’s.1  Let  Af  = l/(NAt)  = 
1/T,  F = NAf.  Define 

gp(t)=  ? g(t  + CT).0<t<T.  (7) 


i.-' 

V • 


h(f  + mF) 4;  < I < 172  . 


(S) 


Theorem  1 (Cooley). -If  h is  the  Fourier  transform  of  g,  then  h is  the  FFT  of  (where  the 
appropriate  sampling  is  made). 


Definition.  Aliasing  is  that  error  introduced  into  the  FFT  by  the  contribution  of  frequencies 
higher  than  those  considered.  Observe  that 

a.  If  g(t)  = 0 for  t > T and  t < 0,  then  g It)  = git).  Ilms,  the  FFT  of  g agrees  with  hp  sampled 
at 


kAf 


n /: 

N/2  * i 


lii  (Ids  case,  the  FFT  evaluated  at  mAf  differs  from  the  Fourier  transform.  It,  evaluated  at  mAf  by 

/ lilmAf  ♦ VF) . 

C/u 

b.  If  ltd*)  - 0 for  i F/2.  then  h(llkAf)  agrees  with  It(kAf)  tor  N/2  ♦ I s k N'2.  In 
this  case,  the  inverse  FFT  sampled  at  nrAt.  0 m ^ N differs  from  g(inAl)  hy 


/ glAt  ♦ V!  T ) . 
k’/o 


h.  H‘.  O hi  tty,  ('.A.  IV.  I.ewtt,  iinj  I1,  l).  Welch,  The  Unite  Ftturict  Ihtwform,  I IFF  7>j«i.  AuJui  Fleettt  umui,..  t_T 
I June  7NSJ5. 
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c.  If  both  the  conditions  for  “a”  and  “b”  are  met,  then  the  FFT  agrees  with  the  Fourier 
transform  at  the  specified  points.  Unfortunately,  the  only  allowable  function  satisfying  these  criteria  is 
g a 0 (*ince  a nonzero  continuous  function  cannot  be  both  band  limited  and  time  limited). 

3.2  Applications  of  Theorem  1 

3.2.1  Choice  of  N and  At  for  Reasonable  Approximations 

Given  g,  choose  S large  enough  so  that  g(t)  is  negligible  for  t > S,  and  choose  N large 
enough  so  that  h(f)  is  negligible  for  |f|  > F/2  where  F = N/S.  Then  let  At  = S/(N  - 1)  . In  this  case,  the 
FFT  is  a reasonable  approximation  to  the  Fourier  transform  sampled  at  these  points. 

3.2.2  Inverse  Transforms 

If  one  is  given  equispaced  samples  of  the  Fourier  transform  and  wishes  to  use  the  FFT  to 
evaluate  the  time-domain  function,  one  should  try  to  duplicate  hp  as  closely  as  possible  by  single  or  double 
aliasing. 

Single  Aliasing. — If  equispaced  samples  are  available  merely  in  the  range  [0,  F/2) , construct 
an  equispaced  array  whose  second  half  is  the  conjugate  of  the  reflection  about  F/2. 

In  FORTRAN  notation,  single  aliasing  is  performed  by 

N1  = N/2+1 

N2  = Nl-2 

F(N1)  = 2*RLAL(F(Nl)) 

DO  1 I = 1,  N2 

1 1 F(N1  + I)  = CONJG(F(Nl-I))  . 

Double  Aliasing.— If  equispaced  samples  are  available  in  the  frequency  range  [0,  F] , replace 
the  second  half  of  the  array  with  the  sum  of  the  original  value  and  the  conjugate  of  the  reflection  about 
F/2.  The  first  half  of  the  array  is  then  replaced  by  the  reflection  of  the  new  second  half.  Double  aliasing 
may  be  achieved  in  FORTRAN  by 

DO  2 I = 1 , N2 

F(N  1 + I)  = F(N  1 + I)  + CON  JG(F(N  1-1)) 

2 F(N  1 -I)  = CONJG(F(N  1 + 1)) . 

For  both  single  and  double  aliasing,  sufficient  information  is  often  not  available  to  com- 
pletely determine  the  best  0-frcquency  value  to  use.  In  the  inverse  transform,  this  can  lead  to  either  drift  of 
the  time-domain  function  or  errors  for  small  values  of  t.  Thus,  often  the  analyst  anchors  the  transform  by 
subtracting  from  the  points  of  the  time  array  the  value  at  0 time.  Doing  so.  however,  may  introduce  a 
relative  displacement  of  the  resulting  time  function,  often  noticeable  at  late  times. 

4.  SUBROUTINL  FLAT 

In  light  of  the  inherent  deficiencies  of  the  FFT  pointed  out  by  theorem  1 and  the  nature  of  the  digiti- 
zation process  employed  by  HDL.  an  alternate  fast  approximation  to  the  Fourier  transform  was  desired. 


I) 


This  transform  was  named  "FLAT”.  It  approximates  the  given  function,  g,  by  a piece-wise  linear  function, 
g2,  and  then  uses  the  Cooley-Tukey  algorithm  to  obtain  a sampling  of  the  Fourier  transform  of  g2  at  the 
same  frequency  points  as  does  the  FFT.*  Approximation  of  continuous  functions  by  line  segments  is  in- 
herently better  for  integration  than  approximation  by  impulses  as  employed  by  the  FFT.  Computer  running 
times  for  FLAT  are  within  25  percent  of  those  for  the  FFT  (approx  200  ms  for  a 2048-point  complex 
array  on  a Control  Data  Corporation  (CDC)  6600  computer). 


4.1  Derivation  of  FLAT 

Consider  a function  g(t)  defined  on  the  interval  [0,  (N  - l)At]  with  g(0)  = O'  extend  g to  the 
point  (NAt,  0)  linearly.  Let  g2(t)  be  the  piece-wise  linear  function  joining  ^ j At,  g()At)J  to  f(j  + l)At, 
g[G+l)At])  , j = 0,  N — 1: 

N-l 

g2(t) = y g2jW  > 

where 

g2j(t)  = cjt + dj  > 
teJjAt,  (j  + l)At] . 

jc  _gQ->-  l)At;~gOAt) 

■ 1 At 


Then  the  Fourier  transform,  G2,  of  g2  is  giver  (for  f 0)  by 


r NAt 

g2(0  =/_  g2(»(|,2"ft|  dt  = l dt 


r (j+l)At 


= i f + 

j=0  ■'jAt 


d« 


_c  J"e(p27rif(j-l)At)  _e(p2ffifjAt)j 


j“0 


<2rri0 


N 


1 I s.e(p2TrifAt) 


(2nif)2  j=0 


PASS. 


*The  Cooley-Tukey  algorithm  used  Is  another  one  of  the  brilliant  creations  ofW.T . Wyatt  of  HDL,  written  In  COM- 
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4.1  Derivation  of  FLAT 

Consider  a function  g(t)  defined  on  the  interval  [0,  (N  - l)At]  with  g(0)  = 0- extend  g to  the 
point  fNAt,  0)  linearly.  Let  g2(t)  be  the  piece-wise  linear  function  joining  (jAt,  g(jAt)J  to  ((j  + l)At, 
g[(j+l)At]j  ,j  =0,N- 1: 

N-l 


g2(t>  = y g2i(t) , 


82jW  = cjtvdj’ 

te[jAt,(j  + l)At}, 


_gQ  + l)At;-gQAt) 
i " At 


Then  the  Fourier  transform,  G2>  of  g2  is  given  (for  f ^ 0)  by 


/■  °°  fNAt 

G2(0=|  g2(t)e(p2nift)  dt  =jg  g2(t)e(p27Tift)  dt 


N-i  r (j+l)At 


1 / 

j=0  Jj 


*i(t),<p2“,dt 


^ (p2?rifO-l)At)  _e(p2nrifjAtjj 


1 V (p27ritAt) 

' - L,  s,e  , 

(2rriQ2  J“0  J 


*The  Cooley-l'ukey  algorithm  used  Is  another  oneof  thebrillLmi  creations  of  W.T.  Wyatt  of  HDL,  written  in  COM- 


where 


U 

fcf 

t; 

£ 

}{ 


\r-1 
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Now  let  Af  = l/(NAt): 


SN  CN-1  ’ 

sj = cj  — Gj-i  » 0<  j < N - 1 . 


. N /p2mjk\ 

Gj(kAo = — — „ I «ri r) 

(27rikAf)'1'  j=0  J 

1 T^1  /p27rijk\ 

= ; X lie\  N ), 

(2rrikAf)2  J=0  f 


where 


' co  CN-1 


1 j = 5 i 


for  j >0. 

Next,  observe  that 


N-l 

I 


j=0 


p27rijk\ 
t.e\  N / 


may  be  evaluated  by  use  of  the  Colley-Tukey  algorithm. 


It  is  of  interest  to  observe  that  if  |g  (t)  - g2  (0  I < e,  0 t <T  , then 

g(t)e(p2,rifl)  dt  — f g2(t)e  (P2mft>  dt 
^ 0 

£ 

% independent  of  f.  This  implies  that,  if  g(t)  is  continuous  and  vanishes  for  t > T,  the  Fourier  transform  of  g2 

| converges  uniformly  to  the  Fourier  transform  of  g as  At  -*  0.  Hence,  the  values  of  FLAT  converge  to  the 

•r  values  of  the  Fou'.er  transform  at  the  sampled  points. 

An  inverse  to  FLAT,  called  “FLIT,”  also  uses  the  Cooley-Tukey  algorithm.  Subroutine  FLIT  produces 
an  equispaced  time  array, 

: ME-'o 

\ with  g(Q)  = 0,  such  that  FLAT  of  the  (linearly  interpolated)  array  is  the  given  equispaced  frequency  array. 


< 


f 


lg(0  - g2(t)  I dt  < eT 
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Indeterminacy  of  the  O-frequcncy  value  may  result  in  a tilt  of  the  time-domain  curve  or  oscillations 
for  large  values  of  t. 

Subroutines  FLAT  and  FLIT  are  resident  in  the  user  library  ANAPAC  on  the  CDC  6600  computer  at 
Fort  Belvoir,  VA. 

The  calling  sequence  for  FLAT  is: 

Call  FLAT  (ARRAY,  N,  DT,  I), 


ARRAY  = name  of  complex  array  containing  data, 

N = number  of  points  in  array  (N  must  be  a power  of  2), 

DT  = At, 

I = value  of  p desired  in  transform. 


The  calling  sequence  for  FLIT  is 

Call  FLIT(ARRAY,  N,  DF,  I), 


ARRAY  = name  of  complex  array  containing  data  (only  the  first  N/2  + 1 points 

need  be  specified);  the  desired  output  is  the  real  part  of  ARRAY, 
N = number  of  points  in  array  (N  must  be  a power  of  2jj, 

DF  = Af,  / 

I = value  of  p desired  in  inverse  transform.  / 


The  desired  output  is  the  real  part  of  ARRAY, 


SUBROUTINES  NUFT  AND  INUFT 


Subroutines  NUFT  and  INUFT  are  used  to  compute  direct  and  inverse  Fourier  transforms. 
The  arguments  in  the  CALL  NUFT  card  are  / 


(1)  Real  array:  independent  input  variable  (time)  / 

(2)  Real  array:  dependent  input  variable  (amplitude)  / 

(3)  Integer:  number  of  points  in  input  arrays  / 

(4)  Integer:  number  of  points  in  output  arrays  / 

(5)  Real  array:  independent  output  variable  (frequency)  / 

(6)  Integer:  indication  of  whether  the  array  in  No.  5 is  given  as  frequency  or  circular  fre- 
quency and  whether  it  remains  that  way  or  changes  (see  comments  for  10M| 

(7)  Complex  array:  dependent  output  variable  (Fourier  transform) 

(8)  Real  array:  storage  for  intermediate  results  to  save  computation  time,  dimensioned  as 

the  input  arrays  / 

(9)  Integer:  sign  of  the  exponential  (+1  or  -1).  / 
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Terms  with  small  absolute  value  in  the  exponential  are  computed  by  use  of  a power  series  expansion, 
to  avoid  loss  of  accuracy  due  to  subtraction  of  almost  equal  numbers. 

The  Fourier  transform  is  computed  assuming  that  the  time-amplitude  trace  is  formed  of  straight  line 
segments. 


It  is  possible  to  use  NUFT  to  perform  an  inverse  Fourier  transform  by  applying  it  separately  to  the 
real  and  imaginary  parts  of  the  given  transform.  This  operation  takes  approximately  twice  as  long  as  the  di- 
rect Fourier  transform.  In  INUFT,  those  operations  are  performed  only  once,  just  as  in  NUFT.  The  call  is 
quite  similar  to  that  of  NUFT,  but  arrays  No.  1 and  2 are  for  the  output,  arrays  No.  5 and  7 are  for  the  in- 
put, and  the  scratch  array  No.  8 has  to  be  complex  and  of  dimension  given  by  the  integer  in  argument 
No.  4.  The  values  of  the  Fourier  transform  are  given  only  for  nonnegative  frequencies,  and  the  output  func- 
tion is  assumed  to  be  real. 

These  transforms  have  to  be  used  when  it  is  not  convenient  to  have  equispaced  input  and  output 
points  in  equal  numbers,  usually  a power  of  2.  It  is  often  important  to  use  frequency  points  with  intervals 
that  increase  with  increasing  frequency,  to  sample  closely  the  low-frequency  values,  and  to  go  to  very  high 
frequencies  at  least  with  a few  points.  The  results  of  inverse  Fourier  transforms  are  improved  by  this 
procedure. 

6.  SUBROUTINES  LINT  AND  CLINT 

Subroutines  LINT  and  CLINT  are  used  to  interpolate  linearly  between  points  in  a given  array.  The 
arguments  in  LINT  are 

(1)  Real  array:  independent  input  variable 

(2)  Real  array:  dependent  input  variable 

(3)  Integer:  number  of  points  in  input  arrays 

(4)  Integer:  number  of  points  in  output  array 

(5)  Real:  increment  for  independent  output  variable 

(6)  Real  array:  dependent  output  variable. 

The  independent  output  variable  starts  with  the  first  value  of  the  input  array.  If  any  points  fall  out- 
side the  range  of  the  input  variable,  a value  of  0 is  given  to  the  output. 

Subroutine  CLINT  differs  from  LINT  only  in  the  sixth  argument,  which  is  a complex  array. 
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7.  MIMIPULSE* 

A model  that  has  been  used  to  represent  analytically  the  response  to  an  electromagnetic  pulse  is  a 
function  that  vanishes  for  t < 0 and  is  defined  by 


, , f nit  - T ,1 

Aexp[— £j(t  — Tj)]  cos  -^r— 


a + 0t  + 6t2  + 7t3, 


aexp[-S2(t-T3)]  cos  [w^t-Tj)] 

+ (a ?2/<o2)  exp  [-|3(t  - T3)]  sin  [w2(t  - T3)] , 


0 < t < T, 


Tj  <t<T2 


T2<t<T3 


T3<t<TE 


t>TE. 


The  values  of  A,  Tj,  T2,  T3,  and  a are  given.  In  addition,  a time,  TE>  is  given  to  define  the  end  of 
the  “recorded”  part  of  the  pulse.  The  values  of  and  £3  are  given  through  constants  Sj , S2,  and  Sjby 


= 3l(T2  ~ Tj) 


and  the  circular  frequencies  are  given  by 


^S2(TE-T3)’ 


3"S3(TE-T3) 


2(T4  - T3) 


2 " TE  " T3 

Finally,  a,  (3,  fi,  and  y are  determined  by  matching  the  values  and  the  derivatives  of  the  function  at  the  ends 
of  the  interval, 


*A  concept  originally  developed  by  Carl  Konschitik  formerly  of  HDL. 
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The  Fourier  transform  of  this  function  is  the  sum  of  the  Fourier  transforms  of  the  components  in  the 
different  intervals: 


Afexpf-koTj^icoTj  + 1)  - l] 


F2(co)  = Aexp^icoTj) 


Tjg/ 


(£j+iw)2  + 


1 -1 


4(T2-T^ 


F3(co)  = exp^-koT3) 


X I “Pl-(*t  + “)(T2  -T.)l  + (S>  +i")l 

26  + 67Tj  67 


■ 3 4 

-ICO  ICO  co 


F4(co)  = aexp(-koT3) 


- exp^-iwT^ 
£2  + ico 


(?2  + ico)2  + (o2 


Y 26  + 67T2  67 

co2  iw3  co4  J 


-+r? 


exp[-(g2-Hco)(TE  -t3)^ 

2x,,2 


($3  + iw)  + ^2  ($  2 + iw)  +w 

X [— (^2  + iw>  cos  wl(TE  ~ T3) + wl  sin  wl(TE  “ T3)] 


^2  exPH?3+H(TE“T3)] 
+ * L-  ■ • ~ ~2 


Wn 


{h+i“Y 


+ C0, 


X [-({3  + iw)  sin  w2{Te  - T3)  - w2  cos  w2(TE  - T3)] 


'] 


where 


7 = - tl  - + : , 

(h~hf  (T!-V2) 

Y = “ 2(T2  - T,j  cxP[”hll2  1 1 *]  ’ 


6 = 


(V1^ 


7(T2  + 2T3) 


and  r;  is  1 or  0,  depending  on  whether  one  assumes  that  the  trace  vanishes  alter  2 or  is  given  by  the  ex- 
pression in  the  last  interval. 


lb 
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Functions  Fj  and  F3  have  to  be  computed  separately  for  cu  = 0,  They  become 

AT, 


F3(0)  = a(T3-T2)  + 


Fl^°)  2 * 

«(Tf-Ti) . (t!-t|)  «K-t1) 


where 


a = -T20  - T26  - T27  , 
j3  = -2T36  - 3T3y  . 

Three  characteristics  of  the  function  follow:  (1)  It  is  continuous  at  the  extremes  of  the  intervals,  but 
it  does  not  necessarily  vanish  at  Tf.  (2)  The  derivative  is  discontinuous  at  Tj , but  continues  at  T2,  where 
the  function  vanishes,  and  at  T3,  where  the  function  reaches  a minimum  given  by  a.  (3)  At  T4,  the  first 
term  in  the  corresponding  expression  for  f(t)  vanishes. 

Subroutines  ANMI  and  ANTRA  compute  the  values  of  f(t)  and  its  Fourier  transform.  Two  real  arrays 
are  passed  to  ANMI,  one  for  the  input  (times)  and  one  for  the  output  (amplitudes),  plus  an  integer  that 
gives  the  number  of  points.  Also,  the  subroutine  reads  a NAMELIST  card  PULSE  that  contains  the  values 
A,  AA,  Tl,  T2,  T3,  T4,  TE,  SI , S2,  S3,  and  S4;  the  meanings  are  obvious,  with  the  possible  exception  of 
AA  = a.  A real  array  for  the  input  (circular  frequencies)  and  a complex  array  for  the  output  (values  of  the 
Fourier  transform)  are  passed  to  ANTRA,  with  an  integer  for  the  number  of  points.  This  subroutine  also 
reads  a NAMELIST  card  with  the  name  PULSE.  It  contains  the  same  information  as  that  for  ANMI,  plus  a 
real  variable  FINIT  that  takes  the  value  0.  for  an  infinite  trace  and  1 . for  a finite  one. 

8.  SMOOTHING 

The  power  spectrum  obtained  from  a digitized  time-amplitude  trace  shows  a strong  oscillating  noise 
component,  especially  on  a logarithmic  scale.  Subroutine  SMUZ  can  be  used  to  present  the  output  in  a 
more  intelligible  form.  The  input  is  an  array  of  a function  given  at  constai  1 intervals.  When  two  maxima  or 
two  minima  occur  closer  than  a prescribed  number  of  points,  the  function  in  between  is  replaced  by  an 
average  value  between  the  straight  lines  joining  the  maxima  and  those  joining  the  minima,  with  a tapered 
beginning  and  end.  The  process  is  repeated  until  there  are  no  changes  in  a complete  pass.  Two  different 
thresholds  can  be  prescribed  for  two  sectors  of  the  function.  The  number  of  passes  through  the  procedure  is 
printed  in  the  output;  a number  of  four  to  six  passes  is  usual.  The  computation  of  the  average  ordinate  is 
performed  by  the  subroutine  AVRG  called  by  SMUZ.  The  parameters  in  CALL  SMUZ  are: 

( 1 1 Real  array:  ordinates  of  the  function  being  smoothed 

(2)  Integer:  number  of  points  in  the  array 

(3)  Integer:  thrc^'old  number  of  points  for  the  first  part  of  the  curve 

(4)  Integer:  threshold  number  ot  points  for  the  second  part  of  the  curve 

(5)  Real:  fraction  of  points  in  the  first  part  of  the  curve. 

The  thresholds  have  to  be  chosen  so  that  the  unwanted  noise  is  eliminated  while  the  significant 
extrema  remain;  a number  to  start  with  might  be  1/100  of  the  total  number  of  points. 
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It  has  been  found  that  the  use  of  this  subroutine  on  the  real  and  imaginary  parts  of  the  Fourier  trans* 
form  that  subsequently  has  to  be  inverted  tends  to  introduce  spurious  oscillations  for  late  times. 

9.  DIGITIZATION  AND  TRANSFORM  ERRORS 

Digitization  is  routinely  performed  at  HDL,  with  a Science  Accessories  Corp.  GP2  digitizer,  which 
allows  for  rear  enlarged  projection  of  oscillograms.  The  tablet  is  20  by  20  in.  and  has  a. definition  of 
0.01  in.  The  output  of  this  machine  is  cards  punched  on  an  026  keypunch.  Normally,  recordings  of  a par- 
ticular event  consist  of  four  or  more  traces  at  differing  sweep  speeds.  Further  processing  of  the  data(e.g., 
time  tying,  sequence  checking)  is  performed  on  a CDC  6600  computer  with  a complete  and  innovative 
software  package  developed  by  T.V.  Noon  of  HDL  (unpublished). 

To  evaluate  the  effects  caused  by  the  digitization  and  transform  process,  numerous  experiments  were 
performed  using  Mimipulse  (sect.  7).  Plots  of  Mimipulse  were  constructed  by  the  same  instructions  to  the 
operator  as  would  be  used  for  oscillograms  processing  (e.g.,  marking  just  the  endpoints  of  apparent  straight 
lines). 


Several  different  results  are  presented  (all  compared  to  the  analytic  transform  of  Mimipulse): 

a.  Transforms  of  equispaced  Mimipulse-to  assess  transform  errors  (fig.  1-8). 

b.  Transforms  applied  to  the  digitized  points-to  determine  the  combined  errors  due  to  digiti- 
zation and  transform  (fig.  9,10). 

c.  Transforms  applied  to  the  array  whose  values  arc  the  actual  values  of  Mimipulse  at  the 
digitized  time  points-to  simulate  sampling  errors  (fig.  11). 

d.  Transforms  applied  to  tire  array  whoso  values  arc  the  truncation  of  the  analytic  values  at 
the  digitized  time  points-to  simulate  quantification  errots  (fig.  12). 

c.  Transforms  applied  to  the  array  whose  values  are  the  truncation  of  the  unalytlc  values  plus 
a random  number  (botween  -10  and  10)  of  points  on  the  digitization  tablct-to  simulate  the  offects  of 
the  digitization-transform  process  (fig.  13,  14). 

One  major  conclusion  that  can  bo  drawn  from  this  set  of  experiments:  At  least  for  oscillograms  of 
traces  similar  to  Mimipulse,  errors  introduced  by  tiie  digitization-transformation  process  arc  minor  in  the 
frequency  range  of  interest  (up  1000  MHz)  compared  to  the  inherent  errors  of  the  typical  data-taking 
apparatus  (where  S percent  errors  are  considered  reasonable). 

Also  plots  are  includ*  1 where  inverse  transforms  applied  to  the  (analytic)  frequency  spectra  of  given 
functions,  two  functions  arc  compared  to  the  original  functions  used: 

a.  Mimipulse  (fig.  15-20} 

b.  Double  exponential,  i.e.,ea*  - (fig.  21-23) . 

One  might  conclude  from  these  results  that  FLIT  performs  quite  credibly.  However,  the  whole  area 
of  inverse  Fourier  transforms  is  fraught  with  danger.  In  the  event  that  the  time  frequency  interval  or  both 


18 


'VW-NSt'Kl* Snr»  'A 


^'t®Rsrr?wewK»5 pt>a»  — . • 


do  not  cover  the  (significant)  domain  of  the  function,  results  from  (1)  the  FFT,  (2)  the  FFT  with  the  value 
at  0 subtracted,  or  (3)  FLIT  can  be  unreliable.  Cases  can  easily  be  constructed  by  truncating  the  time  or  fre- 
quency range  of  an  analytic  function-where  either  (1),  (2),  or  (3)  will  give  the  best  results.  Since  all  three 
converge  (pointwise)  to  the  correct  value  as  the  frequency  range  and  N go  to  infinity,  one  may  increase  the 
number  of  points  in  the  array,  provided  the  data  and  computation  facilities  support  this  increase.  Lack  of 
convergence  is  often  indicated  by  failure  of  FLIT  to  return  to  0 or  of  the  FFT  beginning  significantly  far 
from  0.  Convergence  is  often  indicated  when  the  FFT  and  FLIT  agree.  When  the  data  support  it,  INUFT, 
with  its  ability  to  accept  nonequispaced  frequency  records,  may  be  the  transform  chosen. 


Figure!.  Absolute  values  of  analytic  transform  of  Mimipulsc  (solid  line)  and 
FLAT  of  equispaeed  Mimipulsc  for  256  points  (dotted  line). 


ASJP'-CT<££ 
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Figure  12.  FLAT  of  analytic  values  of  Mimipulse  evaluated  at  digitized  time  points,  quantified  to 
multiples  of  l/l  000  of  maximum  value  and  linearly  interpolated  to  2048  points, 
simulating  sampling  and  quantification  errors  (solid  line)  and  analytic  transform  of 
Mimipulse  (dotted  line). 


FREQUENCY  (x(0«) 

Figure  13,  FLAT  of  analytic  values  of  Mimipulse  evaluated  at  digitized  time  points  quantified  to 
multiples  of  1/1000  of  maximum  value  and  R multiples  of  1/1000  of  maximum  value, 
where  R is  random  number  between  -10  and  10,  simulating  effects  of  digitization  (solid 
line)  and  analytic  transform  of  Mimipulse  (dotted  line). 
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10.  SUBROUTINES  RDTAPE  AND  CSTOUT 

Subroutines  RDTAPE  and  CSTOUT  were  written  by  T.V.  Noon  of  HDL,  and  they  are  used  to  read 
the  digitized  data  for  storage  and  prepare  them  for  input  for  further  analysis.  Subroutine  RDTAPE  reads 
(possibly)  multiple  data  sets  consisting  of  binary  information,  fills  up  the  appropriate  arrays,  and  checks 
for  the  end  of  the  collection  of  data  and  for  irregularities  in  the  data  format. 


by 


Subroutine  RDTAPE  may  be  used  in  two  modes,  only  one  of  which  is  of  interest  here.  It  is  called 


Call  RDTAPE(NT,  X.Y.N, LABEL) 


NT 

X 

Y 

N 

LABEL 


number  of  file  containing  data  (integer) 
name  of  array  to  receive  independent  variable, 
name  of  array  to  receive  dependent  variable, 
name  of  (integer)  variable  to  receive  number  of  points, 
name  of  array  to  receive  label  (eight  words), 


Subroutine  CSTOUT  checks  the  output  of  RDTAPE  for  monotonicity.  If 

X0+  l)  = X(I  + 2)=...  = X(I+n), 

then 

X(I  + 2) X(I  + n),  Y(I  + 2) Y(I  + it) 

arc  removed  from  the  respective  arrays,  the  arrays  are  reordered,  and  Y(l  + 1 ) is  replaced  by 


Y(  I + i)/n. 


II 


then 


X(l  ♦ 2)  < X(1  ♦ I ) , 


X(I  + 2).  Y(l  ♦ 2) 

arc  removed  from  the  respective  arrays,  and  the  arrays  arc  reordered.  The  routine  U called  by 
Cali  CSTOUT(X,Y,N) 
where 

X = name  of  independent  array, 

Y = name  of  dependent  array, 

N = number  of  points -possibly  a new  value  will  be  returned. 
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APPENDIX  A.-USTINGS  OF  PROGRAMS  USED  IN  NUMERICAL  FOURIER  TRANSFORMS 
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SUBROUTINE  FLAT  (YNVQA»N5TAR,DT » JFLAG) 

COMPLEX  X,YNYQA(NSTAR| ,YV1,YV2,XM,TCI ,A 
C FLAT  CALCULATES  THE  TRANSFORM  TO  FREQUENCY  SPACE  OF  THE  COMPLEX  ARRAY 
C GIVEN  IN  YNYQA.  NSTAR  IS  THE  NUMBER  OF  POINTS  IN  THE  ARRAY..  OT  IS 
C DELTA  TIME  FOR  THE  ARRAY.  JFLAG  SHOULD  BE  SET  TO  -1  IF  THE  TRANSFORM 
C IS  TO  BE  EXPRESSED  IN  TERMS  OF  EXP I-2*PI*F) ,+l  OTHERWISE.  THE  TRANSFORM 
C IS  DONE  IN  PLACE. 

YV1  = 0 .0  $ NPT=NSTAR-1  $ SNYQ*0.0 
TPI=8.*ATAN(1.)  1 TP12I* (-1 ./ (TPl*TPI I ) 

THAX=NPT*DT 

Nl=NSTAR/2«-l  » N2=Nl-2 

0F=1./(NSTAR*0TI 

TCI=JFLAG*DF*CMPLX(0.,TPI> 

T=DToNPT 
AxYNYQAINSTAR J 
DO  10  I *2  *NPT 
SNYQ=SNYQ*REAL(YNYQA(  I 1 ) 

10  CONTINUE 

SNYQxSNYQ*  . 5*RE  AUYNYOAI  NSTAR  H 
600  00  HO  1*1, NPT 

YV 1= YNYQA ( I ) 

YV2*YNYQAI 1 ♦ 1 1 
110  YNYOf,  I)*YV2-YV1 
YNYQA(NSTARI*-YV2 
YV  2*0 . 

500  DO  115  1*1, NSTAR 
YV1*YNYUA( I) 

YNYQA!  I I-IYV2--VV1 l/DT 
115  YV2-YV1 

YNYQ A( 1 )«VNYQA( II ♦YVl/OT 
CALL  FF T1D( YNYQA ,NSTAR, JFLAG) 

YNYQ A ( 1 ) *SNYQ*DT 
DO  175  1*2, N1 
X"lTCI*n-l>lo*2 

175  YNYQ A( 1 1 * -YNYQA ( II /X 
00  176  I-I.N2 

176  YNYQA (N 1*1  l«CONJG( VNYQAINl-I I ) 

RETURN 

END 

SUBROUTINE  FLIT  I YNYQA , NSTAR , OF , JFLAG » 

COMPLEX  YNYQA INST AR ) ,U ,S 

C FLIT  CALCULATE  THE  TRANSFORM  TO  TIME  SPACE  OF  THE  COMPLEX  ARRAY  GIVEN 
C IN  YNVQA.  NSTAR  IS  THE  NUMBER  OF  POINTS.  OF  IS  THE  FREQUENCY  INCREMENT 
C JFLAG  SHOULD  BE  SET  TO  *1  IF  THE  FREQUENCY  ARRAY  IS  EXPRESSED  IN  TERMS 
C OF  EXP<-2*P|*F».,  -1  OTHERWISE,  OUTPUT  IS  A COMPLEX -ARRAY ■ YNYQA 
C WHOSE  REAL  PART  IS  THE  DESIRED  TRANSFORM. 

A-RE  AL  < YNYQA 1 1 1 1 
OT  *| ./(DF*NSTAR I 
TP|-8.*ATAN( 1.1 
N1  "N5TAR/2+1 

M2 -Ml -2  ■ - 

FAC*-tTP! *OF I ••37N5TAR 
DO  105  l-l.Nl 

105  YNYQA (I )* YNYQA  1 1 )*FAC*t I -I l**2 

00  110  I -1 ,N2 

YNVQAIN1*! 1-CQNJGI YNYQAINI-1 1 1 


i. 
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110  CONTINUE 

CALL  FFTIDIYNYQA.NSTAR.JFLAG ) 
XL«0  . 

00  75  I=2,NSTAR 
XL=XL+REAL{YNYQAII ) I 
75  CONTINUE 

YNYQAIl 1 =— XL 
U>0.  $ T>0. 

S=0. 

00  50  I =1  ,-NSTAR 
S = 5*YNYQAm 
YNYQ A ( I I »N 
50  N=H>S»0T 

72*0  . 

T=DTo(N5TAR-lI 
X = 0. 

00  325  1 = 1 , NST AR 
325  XaX*REAL(YNYQA(  1 1 1 
N = NS  T AS 

X2=( A/0T-X>*2./IN*IN-1»> 

03  A 35  I =1 »N5TAR 

435  VNYQAm*VNVQA<n«U~ll«X2 
RETURN 
END 
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SUBROUTINE  ANM1(X,V,N) 

DIMENSION  XI N > » Y( N I 

NAHEL 1ST  /PULSE/  A,AA,TltT2,T3,TG,TEfSl»S2,S3,SA 
PI*A.*ATAN(1.) 

READ  PULSE 
PRINT  PULSE 
Wl=PI/CT<*-T3)/2. 

U2»P!*S^/ITE-T3I 
*1=1 «/Sl/(T2-Tl I 
X2=l  ./S2/ITE-T3) 

X3=l ./S3/(TE-T3 ) 

¥Y=-Pl*A/2./(T2-Tl»*EXP  I -X 1 * < T2-T1 ) ) 

GG=-2.«AA/IT3-T2)**3*YV/(T3-T2)*»2 

D&=-AA/(T3-T2I»*2-S&*CT2*2.*T3> 

C1=A/T1 
C21*T2-T1 
C22=PI/2./C21 
CA*«X2/W2 
DO  16  1=1, N 

IFtxm.GT.Tl>  GO  TO  1? 

CONTINUE 

Nl-l-1 

DO  18  1*N1,N 
IFiX(M.GT.T2l  GO  TO  19 
CONTINUE 
N2-I-1 

DO  20  I«N2,N 
[FiXUI.GT.T3l  GO  TO  21 

continue 

N3-I-1 

DO  22  I •N3 ,N 
IFtxm.GT.TEI  GO  TO  23 
CONTINUE 

NA«N 

GO  TO  2G 

NG«I -[ 

03  30  t *1 ,N1 
CC21«A*EXPUI*T1> 



DO  AO  I*Nt,N2 
Tl-Xill 

till «C C2I#E XP ( -XI *T1I*C0SIC 22*1 TI -Till 
K2-N2M 

8G«i-2.*0G-3.*GG*T3»*T3 

»G«i-BG~t06*G0*T2l*Y2»*f2  •••—•  - 

00  SO  I *N2»N3 

Ti-xm 

Till *AG*(8G* (OG*GG*T| |*T 1 1 *t ! 

NS-NSO 

CC«1*AA*EXI’<X2«T3> 

CCAG»«!*T3  

CC4*.«AA*CAG*EXPiX3*T3> 

CU3«U2*T3 
00  GO  I «N3 »N 
TI »X  1 1 1 

VIII»CCGl*eXPI-X2*TII*C0SiNl*Tl-CCGG»*LCA5nXPl-X3*Tn 
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1 *S1N(W2*TI-CC48I 
IF (N4.EQ.N ) RETURN 
DO  70  I*N4,N 

70  vm«o. 

RETURN 

END 

SUBROUTINE  ANTRA<OH,FT»NPTS> 

COMPLEX  FT(NPTSI,FTI,CWT1,CWT3,CWT4,CH2,CH3,EVE 
COMPLEX  CWT,CWT2,FTI,FT2,FT3,FT4#CM1 
DIMENSION  OMINPTS) 

NAMELIST  /PULSE/  A.AAtT 1 ,T2 ,T 3,TA,TE , SI »S2 , S3.SA.F INI T 

PI*<..»ATANCl.i 

EVE=CMPLX(0. >1 .) 

READ  PJi.SE 
PRINT  PULSE 
Ml=PI/(T4-T3I/2. 

M2=Pt*SA/(TE-T3» 

X 1 =1 ./S1/IT2-T1 ) 

X2-1 »/S2/ I TE-T3 ) 

X3=l ./S3/ITE-T3I 

Y*-PI«A/2./(T2-Tl»*EXP  I -XI *02-7111 

0&"-2.*AA/(T3-T2I**3*Y/(T3-T21*°2 

0G*-AA/IT3-T2I**2-G&*IT2*2.*T 31 

BG«{-2.*0&-3.*GG*T3»eT3 

•&“( -BG-(0G*GG*T2I*T2J*T2 

CO  Is  T3-T2 

Cl *A/Tl 

C21«T2-TI 

C22-PI/2./C21 

C?3«  C22**2 

C2A*EXP(-XI*C2U 

C 3 * * 2 . *0&*6 .°GG*T  3 

C32-2.*0G*6.*G&*'? 

C33*f>.*G& 

C<>  l»TE-T3 

t<x>aWI**2 

CA3«M2*®2 

C*A»X2/W2 

CNS*COS(Ul*C«l I 

C«.6«SIN(MI*£AI)*Ml 

CW*StN(M2*C<>l  > 

C<.8«C3.S(W2*C<in*tl2 
C<.9«E*Pt-X2*CU  I 
C<.I0*EXPC-X3*CAU*CA«. 

DO  t.QO  i-l.HPTS 

MT  »ti  *T  I 
MS«M**2 
MO«MS*M 
MF«tfQ*M 
MH«M*C2l 
MT  2«U*T 2 
MT  J»M*T3 
WT<.*M*C<tl 

CMT»CMPLK<CaS«MT),-SIMI«T! > 
CMU-CMPLX4C0S«um,-SIN(MTn) 
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C*T2*C1PIX(C0S(RT2>,-SIN{WT2)  ) 
CHT3*CHPLX(CDSJ«T3>,-SlNIWT3n 
C»»T4*CXPIX(COSCWT<»),-SINOIT41  > 

Cm*CNPLX(Xl.M) 

CH2=CNPLX(X2,« ) 

CK3*CHPLX(X3.M) 

IF (M .NE .0.0)  GO  TO  92 
FT  I «.5*»<>T1 

FT l=FTI»(AG»9G*(T3*T2)*.5+DG*IT3**2*T3*T2*T2**2)/3.»GG* 
1 <T3**3*T3*o2*T2*T3*T2««2*T2**3»*.25»K0l 

GO  TO  94 

92  FT  I -CUV  “CHPLX  < 1 . tUT  >-l 

FT  I =C 1 *FT I /M S 

FT1*AA«EYE/R-EYE*C31/WQ-C33/MF 
FT  2= Y/rfS-EYE*C32/WQ-C33/«F 
FU«FTl*CWT3*FTl-CMT2*fT2 
94  FTKNl»r.Nl*C23 

FT2«C24«CWT1  *C22*CH1 

FT1«A*CWT  /FT1 

FTI*FTI*FT1«FT2 

F T 1 • AA*CWT3 

FT2--CH2*C45*C46 

FT2  = C49*CMT4*FT2*F  Ui T*CM2 

FT3«-CH3*C4?-C48 

Fr3«C410*CWT4*FT3*FJNi <*C44*«2 

FT<,-CH2*CH2»C42 

FT2»FT2/FT4 

FT4«CH3*CH3*C43 

FT3»m/FT<. 

100  Frm-FT|*FT|*«fT2#FO» 

return 

ENQ 
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SJ8R0UT1NE  NJFT  (X,  Y,N,»U,OM,ION,FT  ,DX  .JFLAG) 

C THIS  SUBROUTINE  CALCULATES  THE  FOURIER  TRANSFORM  OF  THE 
C FUNCTION  GIVEN  8V  STRAIGHT  LINES  JOINING  THE  POINTS  GIVEN  BY  THE 
C ARRAYS  X,Y.  THE  NUMBER  OF  INPUT  POINTS  IS  N,  THE  FREQUENCY  ARRAY  OH  IS 
C PROVIDED  AND  HAS  DIMENSION  NU.  DX  IS  A REAL  ARRAY  OF  DIMENSION  N 
C THAT  IS  USED  FOR  SCRATCH.  THE  FOURIER  TRANSFOR  IS  GIVEN  IN  THE 
C COHPLEX  ARRAY  FT. 

C IF  I OH  ~L  , INPUT  AND  OUTPUT  QM  ARE  CIRCULAR  FREQUENCIES 
C IF  IOH=2, INPUT  ARRAY  IS  CIRCULAR  FREQUENCY,  OUTPUT  IS  FREQUENCY 
C IF  10M=3,  INPUT  ARRAY  IS  FREQUENCY,  OUTPUT  IS  CIRCULAR  FREQUENCY 
C IF  I ON  - A , INPUT  AND  OUTPUT  QM  ARE  FREQUENCIES 
C JFLAG-  *1  IF  THE  FOURIER  TRANSFORM  HAS  A FACTOR  £XP(*IWT) 

C JFLAG=-1  I F "THE  FOURIER  TRANSFORM  HAS  A FACTOR  EXP(-iMT) 

DIMENSION  X(N),Y<N),OX(N>,OHINU) 

COMPLEX  EYE , EX1 ,EX2,S1 ,S2,FT (NU) 

TJPI=8.®ATAN(1.) 

TP !N*l ./TOPI 
S-O. 

EYE* ID.  ,1  . ) 

IF  HOH.LE  .2)  GO  TQ  25 
03  2D  I - 1 ,NJ 
20  OH  ( I I =v)M(  I )»TQP  I 

25  NPaN-l 

03  30  1=1, NP 
33  DX(1 )*X(!«ll-X<!  > 

N3*I 

IF  (OHUI.NE.O.)  GO  TO  <.2 
NO  *2 

33  SO  1*1, NP 

no  s«s*mt*n»v(m«Qxm 

FT  m*;MPLX(S/2,0.0» 

A2  03  AS  1*1 ,NP 

AS  QX(I  I « t Yl  1*1  J-'V(I  ll/OXIi  1 

03  TJ  I *NO»NU 
h*(|H(  1 ) 

#C  *rf S*A 

32*0*3,0.0) 

IF  LAG* 3 
« T * <i  • t I I ) 

EXt*CHPLX(CaS(U(),S|N(«TM 

SI«YUI*tXt 

03  50  J*2.N 

IF  IIFlAG.EQ.il  03  T 3 A? 

*?«*(JI 

Kl«8lJ-ll 

IF (lu*<2t.GT . S.Of -2)  GO  TO  Afc 

»P*X1»L2 

*Q«XP»X1»X2**2 

«*nm**?**i 

OY«YIJ»-YI J-| | 

XRE* I-3.5«mS*XP*«F***/2A.I*D» 

S2*S2*CHPLX(KRE ,XIMI 
GO  f 3 SO 

Sfe  ESI«CHPLX(C3S(tf*Xl),SlN(w**n  1 
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I F Lfi  G = 1 

57  UT  =H  *X  ( J > 

EX2=  CSPLX ( C3S  «WT  » .SINIWTJ) 

S2=$2*iex2-exi)*oxij-i) 

EX1=EX2 

50  C3NTINJE 

IEUFLAG.E0.1I  GO  TO  58 
EXl  = CHPlX»C3S(tm,SIN(WTn 

58  SI  =S  1-Y  IN  I *E  X 1 

70  FTU  >=(EYE*S1*S2/H>/W 

IF  I13N.EQ.1.3R.10M.EQ.3)  GO  TO  35 
00  80  1=1, NU 
80  01 (I )«3M(l  loTPJN 

65  IF  IJFLAG.EO.H  RETURN 

03  90  1=1, NJ 

90  FT  (I  )»C0NJGtFTim 

RETJRN 
END 

SUBROUTINE  I NUFT ( X , Y ,N ,NU ,0N, I ON »f T ,C SCR , JFL AG > 

C THIS  SUBROUTINE  CALCULATES  THE  INVERSE  FOURIER  TRANSFORM  OF  THE  FUNCTION 
C GIVEN  BY  STRAIGHT  LINES  JOINING  THE  POINTS  GIVEN  BY  THE  COMPLEX  ARRAY  FT 
C AS  A FUNCTION  UF  THE  REAL  ARRAY  OH.  THE  NUMBER  OF  INPUT  POINTS  IS  NU. 

C THE  INVERSE  F.T.  IS  ASSUMED  TO  BE  A REAL  FUNC1I0N  AND  GOES  IN  THE 
C ARRAY  V,  ANO  THE  REAL  ARRAY  X IS  GIVEN  ANO  HAS  DIMENSION  N. 

CSCR  IS  A COMPLEX  ARRAY  OF  OIMENSION  NU  AND  IS  USED  FOR  SCRATCH 
C IF  I 0M*1 , INPUT  AND  OUTPUT  OH  ARE  CIRCULAR  FREQUENCIES 
C IF  I0M«2, INPUT  ARRAY  IS  CIRCULAR  FREQUENCY,  OUTPUT  IS  FREQUENCY 
C IF  I ON  *3 , INPUT  ARRAY  IS  FREQUENCY,  OUTPUT  IS  CIRCULAR  FREQUENCY 
C IF  1C1H-5,  INPUT  ANO  OUTPUT  OH  ARE  FREQUENCIES 

C JFLAG«  *|  IF  THE  INVERSE  FOURIER  TRANSFORM  MAS  A FACTOR  EXP(UMT) 

C JFLAG* -I  IF  THE  INVERSE  FOURIER  TRANSFORM  HAS  A FACTOR  EKPI-IWT) 
OIMENSION  XINI .YINI .CMINUI 

COMPLEX  EYE«EXl,£X2»5ltS2»PTI*iUt,£SCR{NUt»S»0V»KHE*XiM  • 
TJPt«B.»ATANU.> 

TP  IS  * * ./TOPI 
P l N»  2 «*TP | N 
S»0. 

E YE*CMPLX 1 0. , 1 .) 

IF iJFlAG.eu.l | GO  TO  15 
03  10  1«I,NJ 
10  FTm-CONJGIFTKI) 

15  tFII3H.LE.2I  GO  TO  25 

03  20  !"! ,NU 
23  ON  < I > »3Mi I |*TOP | 

25  NP-NU-1 

NO®  I 

IFI  xm.NE.O.)  GO  TO  52 
N3  *2 

03  50  l-I.NP 

50  S«S»|FM  1*1)  *FT  I U l«IOM<l*n-OMlin 

VU»«RFAMS»*TP|N 
52  03  55  I-l.NP 

55  CSCRt  SI»IF  f U*n-FT<m/IOHU*l»-OMI  II) 

03  70  t»NO,N 
Hill 

VS*V *v 
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46 
4 7 


50 

48 

70 


80 

85 


90 


«C=W50W 
*iF=WC*K 
S2=(0. 0,0.0) 

IFLAG=0 
tf  T =W  ®0M ( I ) 

EXl=CrtPLX(COSCHT>,SIMCWT>) 

S1=FT(1)*EX1 

DO  50  J*2,NU 

IF  (IFLAG.EQ.l)  GO  .0  47 
X2*=0M  ( J ) 

XI =0M { J-l ) 

IF<<W*X2) .GT.5.0E-2)  GO  TO  46 

XP=Xl«-X2 

X0=XP*XH-X2«*2 

XR=XQ«‘X?.+X2->*3 

OV  =F  T ( J )-FT( J-l ) 

XRC='-0.5*W5*XP+WF*XR/24„)*DY 

xiH*(w-wc*xa/b.)*OY*evE 

S2=$2+XRE+XIM 
GO  TO  50 

EX!=CMPLX(C0S«K*XI),SIN{«*X1) ) 

IFLAG=1 

KT=ri'>0'l(  J) 

EX2  = CMPLX(CQS(WT)  ,51 WtWT  > > 

S2  = S2+(cX2-EXl  )*CSCR(J-1 ) 

EX l1 EX  2 
CONTINUE 

IF (IFLAG.EO.l)  GO  TO  48 
EX1=CRPLX(C05(WT)  ,5  I N ( WT ) ) 
S1=S1-FT(NU)*EX1 
S =(EVE*Sl+52/*n/W 
¥(  I ) -REAL ( S) ®PIN 
CONTINUE 

IF (I0H.EQ.1.0R.I0M.EQ.3)  GO  TO  85 
DO  80  1 = 1, NJ 

QH<n=oNinnpiN 
IF (UFLAG.EQ.l ) RETURN 
DO  90  1 = 1, NU 
FT (I ) =CON JG(FT ( I ) ) 

RETURN 

END 
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SUBROUTINE  CLINT!XF,YF»JF,N5TAR »XNYQ ,YNY0 ) 
DIMENSION  XFIJ?),YFIJF> 

COMPLEX  YNYO(NSTAR) 
bl 

X1=XF ( 1 ) 

YNYQ 1 1 ) =YF (1 ) 

CO  20  L=2,NSTAR 

x=xi*(l-i)*xnvq 
10  IFIX.LE.XFIIJ J CO  TO  20 

1=1*1 

IF(I.CT.JF)  GO  TO  30 
OENOM=XF(I)“XelI-l) 

Cl  = ( YF  ( I ) -YF  1 1 — 1 1 J/DENL.M 

C2=(YF(I— l)*XF(I)-YF(I>*XF(I-i))/  3EN°M 

GO  TO  10 

20  YNYQ{L)=C1*X+C2 

GO  TO  100 

30  DO  CO  J=L,NSTAR 

CO  YNYQ(J)=0. 

100  RETURN 

END 

SUBROU.  INF  LI*)T(XFfYF»JFviNSTAR#XNYB*YNY9> 
DIMENSION  XF( JF) * YF ( JF ) 

DIMENSION  YNYOCNSTARI 
1 = 1 

X1=XF(1) 

YNYQI1 ) *YFd ) 

DO  20  L”2 »NST AR 
X=X1*(L-1 J*XNYQ 
10  tFU.LE.XFdH  GO  TO  20 

1=1*1 

IFd.6T.JF>  GO  TO  30 
DENOM*XFin-XFCI-l> 

Cl=( YFC l>-YF< I— 1) >/OENCHf 
C2=(YF(I-1)*XF(I) -YF (I)*XF{I-1) ) /DENOM 
GO  TO  10 

20  YNYQ (L ) =C 1*X*C2 

GO  TO  100 

30  OQ  CO  J*L  rNSTAR 

CO  YNYOIJl-O. 

100  RETURN 

END 
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SJBROJTINE  RDTAPE  (NT  , X ,Y  , N, LABEL  ) 

DIMENSION  X ( 1 } ,Y(1),LABEL(8) 

READ  (MT)  (LABELS  I I »I  = leB) 

IF  < EOFINT  ) >40,10 
10  READ (MT)N 

I F { E UF ( NT ) ) 50  1 2D 
20  READ{NT)CX(I),Y(I),I=ltN> 

:IF  {EOFINT  J 150,30 
30  RETJRN 
<10  PRINT  41 

41  FDRMAT(5X,*E0F  ENCOJNTERED  PROPERLY*) 

GO  TO  50 

50  PRINT  51, LABEL 

51  FORMAT (5K,*50F  ENCOUNTERED  IMPROPERLY  DURING  *»8A10,*  EXIT  CALLED* 
* ) 

60  CALL  EXIT 
END 


SUBROUTINE  CSTOUT (XF,YF, JF  ) 

cccccccccccccccccccccccccccccccccccccccecccccccecccccccccccccccccccccccccccccccc 


c c 

c CHECK  TIME  ORDERING  OF  THE  POINTS  AND  CAST  OUT  THOSE  POINTS  NOT  IN  THE  C 
C PROPER  TIME  uRDER  C 
C IT  ALSO  AVERAGES  THOSE  POINTS  WHICH  HAVE  THE  SAME  X VALUES  C 
C C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
DIMENSION  XF ( JF ) , YF ( JF ) 

J3  AO  >=0 
J=1 

L SUM  = 1 
K = 2 

XS AV  = XF ( 1 ) 

YSUM  = YM1  ) 

80  IF  CXFUS-XSAV)90,10Q,110 

90  K = K+1 

JB AD  = JB  AD  + 1 
IF  (K -J F ) 80  ,80 , 1 10 
100  L SUN  =L  SUM-fl 

YSUM  = Y SUM  + YF { K } 

K = K+  1 

IF (K-JF )80  ,80  » 30 
110  YF (J }*YSUM/LSUH 
XF  ( J ) =X iAV 
IF(K-JF)120, 130,30 
120  XS  AV  = XF ( K ) 

VSUM>=YF<K) 

L SUM  = 1 
K*K»1 
J=  J+l 
GO  TO  80 
130  J*J*1 

YF(J)=YP<K) 

XF ( J I *XF (K  ) 

30  PRINT  1 , J6AD , JF 

1 FORMAT!/, IX, 15,*  POINTS  OUT  OF  A TOTAL  OF  *,I5, 

1 * DID  NOT  SATISFY  THE  TINE  ORDER  CRITERION*) 
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J3  AO  =0 
1 = 1 

IFtXF(I).GE.O.O)  GO  TO  202 
1 = 1+1 

JB AD  = J8 AD+1 
IHI-JF)  201*205,205 
IF(i.E3.1)  GO  TO  206 
1 = 1-1 
JF  =JF-I 
DO  203  K=1 » JF 
XF IK ) =XF (K+I ) 

YF IK  )=YF(K+I ) 

1 = 1 + 1 

PRINT  210,1 
GO  TO  209 
JF  =1 

PRINT  211 
GO  TO  209 
PRINT  212 
PRINT  213, Jf 
RETURN 

FORMAT! IX, 110  »*  PQIMTS  HAD  NEGATIVE  TINES*) 
FORM ATI IX  **  ALL  POINTS  HAD  NEGATIVE  TINES*) 
FORNATUX,*  NO  POINTS  HAD  NEGATIVE  TIMES*) 
FORM AT(1X, 110,*  POINTS  WILL  BE  USED*) 

END 


kf 


APPENDIX  A 


SUBROUTINE  EFT CA,NP3K ,NST AR ,XNNra , IS  I GN > 

C A = COMPLEX  INPUT  AND  OUTPUT  ARRAY. 

C NPOW  = ANYTHING.  IT  IS  NOT  USED  ANO  IT  IS  KEPT  TO  MAKE  THE  CALL 
C COMPATIBLE  KITH  A PREVIOUS  VERSION  OF  FFT . 

C NSTAR  = NUMBER  OF  POINTS  IN  ARRAY. 

C XNYO  = INCREMENT  IN  THE  INDEPENDENT  VARIABLE. 

C ISIGN  = SIGN  IN  THE  EXPONENTIAL  OF  THE  FOURIER  TRANSFORM. 
COMPLEX  A(NSTAR) 

CALL  FFTIDtA, NSTAR, ISIGN) 

DO  10  1=1, NSTAR 
A ( I ) =A ( I ) *XNYQ 
10  CONTINUE 
RETURN  $ END 
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$5  SB  3 B34-A0 
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RXb 

X1/X5 

SAb 

THETA 

SA1 

THETA 

RJ 

SIN. 

FXI 

Xb»X6 

NX7 

XI 

SA7 

WST 

SXO 

BO 

FXb 
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Xb 
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M 

SB1 

XS 

SA3 

N 

S&b 

X3 

SAO 
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SUBROUTINE  SNUZ ( X , N ,M 1 ,H2 ,FRAC » 

C THIS  SUBROUTINE  AVERAGES  OUT  A FUNCTION  BETWEEN  PEAKS. 

C X IS  AN  ARRAY  OF  N EQUISPACED  ORDINATES. 

C HI  OR  M2  ARE  INTEGERS  THAT  GIVE  A MAXIMUM  NUMBER  OF  POINTS 
C BETWEEN  TWO  MAXIMA  OR  TWO  MINIMA  FOR  THE  AVERAGING  ACTION  TO  TAKE  PLACE. 
C FRAC  IS  THE  FRACTION  OF  THE  POINTS  FOR  WHICH  Ml  IS  USED* 

C THE  REMAINDER  USES  M2. 

C THE  PROCEDURE  IS  REPEATED  UNTIL  NO  CHANGES  ARE  INTRODUCED  ANTWHEREt 
C THE  PRINTOUT  STATES  THE  NUMBER  OF  PASSES  OF  SHUT  THAT  WERE  REQUIRED. 
DIMENSION  XINl 
NT  =0 

10  KFLAG-0 
NT=NT*1 
1FLAG=0 

Nl*3  $ K2»N*FRAC  S H-Hl 
IF4FRAC.NE.0.)  GO  TO  15 
N2*N  i H = N2 

is  xu*xm 

XI =X (2! 

IFtKI.GT.Xll>  GO  TO  20 
LX  =-l 
X2NAX*X I 1 
K2HAX-1 
K2HIN— H 
GO  TO  30 
20  LX “I 

X2HIN-XI1 

K2MIN-1  

K2HAX— N 
30  XU-XI 

AO  00  1000  I»N!,N2 

XI-XIII 

IFfXI.GT.Xm  GO  TO  GO 

IF<LX,EQ.~1»  GO  TO  900  * -•  

LX  — I 

IF  ( I -K2MAX.LE.»M ) GO  TQ  100 
X2NAX«XII 
K2NAX*I -l 
GO  TO  990 

60  IF  (LX. EO. II  00  TO  900  - — — 

LX«I 

IF (I -K2MIN.LE.MI  GO  TO  200 
XZNIN*X  I S 

K2NIN-I-I  * 

GO  TO  990 

100  K1HAX-K2MAX  * — 

XiHAX«X2HAX 

*t2MAX«l“l  - * 

X2NAX-KU 

IFCtftAO.EO.l)  00  10  ISO  - 

KFLA5»*1 

*FUG*i  ■ ■ - ■—  - - — — 

KIN! N*K2MIN-N 

tFIKIMIN.LE.OIKIMN-i  ■ -*■ — — — 

XlHfN-K(KlMIN) 

If  tM  IMAK.EQ.  I > 60  fO  ISO  - — 

CALL  AVRfc<XlMm,KlHIW,X2MI*,«2MIN,XlMI«|,KlHlN#XtMAX#KtIUI,X, 
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1 K1HIN*1 .K1HAX) 

ISO  CALL  AVR&tXlNAX,KlHAX,X2HAX,K2NAX,XlKlM*KllllNtX2HlN,K2HIN,Xt 
1 KlMAXn,K2HlN» 

GO  TO  990 
200  KiMI N*K2N  IN 
X1HIN-X2HIN 
K2HIN=I-1 
X2NIN«XIl 

IFdFLAG.EQ.lt  GO  TO  250 

KFL*G«1  

IFLAG=1 
K1HAK*K2HAX-N 
IFtKIHAX.LE.OI  K1NAX-1 
XlNAX*XtKlNAXl 
IFIKININ.EQ.1)  GO  TO  250 

CALL  A VRGiXIHAX .K1NAX .X2HAX .K2MAX .X1HAX  .K1HAX .XIHIN.KIH IN ,X  , 

1 K1HAX*1 .K1HINI 

250  CALL  AVRG(X1HAX,K1NAX,X2HAX,K2HAX,X1N1N,K1NIN,X2MIN,K2N1N»X, 

1 K1HINK1.K2HAX) 

GO  TO  990 

900  IFdFLAG.EO.OI  GO  TO  990 
IF (LX.EQ.l  I GO  TO  950 
IF  d -K2NIN .LT .N ) GO  TO  990 

CALL  AVR&IXINAX,KINAX,X2NAX,K2HAX,X2NIB,K2HIN,X1,I,X. 

I K2HINM.K2NAXI 

CALL  AVRGtX2NIN,K2NlN,Xl,I,  X2HAX.K2HAX.Xl ,1 ,X,K2MAX*1 ,1-11 

GO  TO  980 

950  IF  d -K2HAX .LT .H ) GO  TO  990 

CALL  A9RGtX2NAX,K2HAX»XI  .1  »X1H1N,K1NIN,X2NIN,K2N1N.X, 

1 K2HAX*l  .K2NINI 

CALL  AVRG(X2KAX,K2HAX,Xl,I,X2NINtK2HIN»XI,l ,X,K2NIN*1 ,1-1 1 
9B0  IFLAG«0 
990  XU«X1 
1000  CONTINUE 

IFtH2.E0.NI  GO  TO  1010 
Nl-N2*l  A N2-N  » H-H2 
GO  TO  AO 

1010  IFdFLAG.EO.OI  GO  TO  1060 
IFtLX.EQ.il  GO  TO  1050 

CALL  AVRGt  X1HAX .K1NAX.X2MAX .H2HAX.X1NIN.K IHIN,X2H|N,K2MIN«X, ■ 

1 K IN  IN* I ,N  1 
CO  TO  1060 

1050  CALL  AVRGtXIHAX,KlNAX,X2NAX.K2NAX,XlHIN'K|HlN,X2NIN,K2NINvX, 

1 KINAXt’  ,NI 

1060  IFtKFLAG.EQ.il  GO  TO  10 

PRINT  1,  N1.H2,  FRAC.NT  

1 FORHATtlX,  >Nl«**,I3,"  N2«-,13,-  FRAC—.F5  .2  , 

1 - NUMBER  OF  PASSES  OF  SHU2  IS-.IAt 
RETURN 
ENO 

SUBROUTINE  AVRGtX) ,KI . X2 ,H2,X 3,K3,X4,K6,X,MI ,KF I 

C THIS  SUBROUTINE  SUBSTITUTES  POINTS  IN  THE  ARRAY  X BETWEEN  Kt  AND  KF 
C BY  THE  AVERAGE  BETWEEN  THE  TWO  STRAIGHT  LINES  DEF1NE0  BY  11.12 
C AND  X3.XA  AT  .POINTS  Kl,  K2,  K3,  HA  RESPECTIVELY 
01HENSIQN  Xtll 
A1-.WIK2HUI 


